knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
# library(sjPlot)
library(multcompView)
library(modelbased)
library(hrbrthemes)
library(patchwork)
library(here)
library(readxl)
library(lmerTest)
library(blme)

# oyster length/weight data
sizdat <- read.csv(here('data/raw', 'field_data.csv'), stringsAsFactors = F) %>% 
  select(Exposure, Organism, Full_ID, Station, Group, Individual, Length_mm, Final_Weight) %>% 
  mutate(
    Station = factor(Station, levels = c('OUT', 'EDG', 'MID')), 
    Exposure = factor(Exposure, levels = c('Pre', 'Post'))
  ) %>% 
  rename(
    Weight = Final_Weight
  ) %>% 
  filter(!Organism %in% 'Mussels') %>% 
  gather('var', 'val', Length_mm, Weight)

# oyster dissolution data first sheet
dissdat1 <- read.csv(here('data/raw', 'Scoring data_RM_NB (1).csv')) %>% 
  rename(Exposure = Station) %>% 
  mutate(
    Exposure = factor(Exposure, levels = c('OUT', 'EDG', 'MID'), labels = c('Out', 'Edge', 'Mid'))
  )

# oyster dissolution data second sheet
dissdat2 <- read_excel(here('data/raw/Kelp Forest SEM ScoringNB.xlsx')) %>% 
  rename(Exposure = Station) %>% 
  mutate(
    Exposure = factor(Exposure, levels = c('OUT', 'EDG', 'MID'), labels = c('Out', 'Edge', 'Mid'))
  ) %>% 
  rename(`Average.Score` = `Avg Score`) %>% 
  filter(!is.na(`Average.Score`))

# combine the two files
# make a unique replicate columns for two replicates at each location/species
dissdat <- bind_rows(dissdat1, dissdat2) %>% 
  unite('newgrp', Exposure, Group, remove = F) %>% 
  mutate(newgrp = factor(newgrp, labels = 1:length(unique(newgrp)))) %>% 
  select(-Rep)

# size models
sizmod <- sizdat %>%
  mutate(
    Station = case_when(
      Station %in% c('EDG', 'OUT') ~ 'Edge/Out', 
      T ~ 'Mid'
    )
  ) %>% 
  group_by(Organism, var) %>% 
  nest() %>% 
  mutate(
    lmmod = map(data, ~aov(val ~ Exposure * Station, data = .x)),
    phmod = map(lmmod, TukeyHSD),
    ltmod = map(phmod, function(x){
      modhsd <- x %>% 
        .$`Exposure:Station` %>% 
        data.frame
      pvals <- modhsd$p.adj
      names(pvals) <- rownames(modhsd)
      lets <- multcompLetters(pvals)
      
      return(lets)
      
    }),
    plors = pmap(list(Organism, var, lmmod), function(Organism, var, lmmod){
      
      sjPlot::plot_model(lmmod, type = 'int', mdrt.values = 'meansd') + ggtitle(paste(Organism, var, sep = ', '))
      
    })
  )

# dissolution models
dissmod <- dissdat %>% 
  group_by(Organism) %>% 
  nest %>% 
  mutate(
    lmmod = pmap(list(data), function(data){

      # out <- lm(Average.Score ~ Exposure, data = data)
      tomod <- data %>% 
        mutate(
          newgrp = fct_drop(newgrp),
          newgrp = as.numeric(newgrp)
        )
      
      # out <- blmer(Average.Score ~ Exposure + (1|newgrp:Exposure), data = tomod)
      # out <- lmer(Average.Score ~ Exposure + (1|newgrp), data = tomod, REML = F) # this is the same as above if all newgrp were unique
      # out <- lmer(Average.Score ~ Exposure + (1|Group:Exposure), data = tomod, REML = F) # this is also equivalent
      out <- lm(Average.Score ~ Exposure, data = tomod)
      
      return(out)
      
    }),
    anomod = map(lmmod, function(x){
    
      out <- anova(x)
      
      return(out)
      
    }), 
    plomod = pmap(list(Organism, lmmod), function(Organism, lmmod){

      # estimates
      mnsval <- estimate_means(lmmod, 'Exposure')
      cnsval <- estimate_contrasts(lmmod, 'Exposure') %>%
        mutate(
          sig = ifelse(p < 0.05, 'sig', 'notsig'),
          sig = factor(sig,levels = c('notsig', 'sig'), labels = c(' not significant', 'significant'))
        ) %>%
        unite('Contrast', Level1, Level2, sep = '-')

      # sample size
      n <- lmmod$model %>% nrow
      # n <- lmmod@frame %>% nrow
      
      # labels
      subttl <- paste0(Organism, ' oyster')
      
      ttl1 <- NULL
      ttl2 <- NULL
      ylab1 <- 'Estimated means (+/- 95% CI)'
      ylab2 <- 'Estimated differences (+/- 95% CI)'
      captns <- NULL#paste0('Significance is where CI does not include zero, alpha = 0.05, total n = ', n)
      if(Organism == 'Pacific'){
        ttl1 <- '(a) Treatment estimates'
        ttl2 <- '(b) Treatment differences'
        captns <- NULL
        # ylab1 <- NULL
        # ylab2 <- NULL
      }
      
      # mean esimate plots
      p1 <- ggplot(mnsval, aes(x = Exposure, y = Mean)) + 
        geom_point(size = 3) + 
        geom_errorbar(aes(ymin = CI_low, ymax = CI_high), colour = 'black', linewidth = 1) + 
        labs(x = NULL, y = ylab1, title =  ttl1, subtitle = subttl) + 
        theme_ipsum() + 
        theme(
          panel.grid.major.y = element_blank(),
          panel.grid.minor.y = element_blank()
        ) + 
        coord_flip()
      
      # contrast plots
      p2 <- ggplot(cnsval, aes(x = Contrast, y = Difference, colour = sig)) +
        geom_point(aes(colour = sig), size = 3) +
        geom_errorbar(aes(ymin = CI_low, ymax = CI_high, colour = sig), linewidth = 1) +
        labs(x = NULL, y = ylab2, title = ttl2, subtitle = subttl,
             caption = captns) +
        theme_ipsum() +
        theme(
          legend.title = element_blank(),
          legend.position = 'bottom',
          panel.grid.major.y = element_blank(),
          panel.grid.minor.y = element_blank()
        ) +
        geom_hline(yintercept = 0, linetype = 'dotted', size = 1) +
        scale_colour_manual(drop = F, values = c('black', 'tomato1')) +
        coord_flip()

      list(p1, p2)
      
    })
  )

Length and weight

plofun <- function(x){
  
  psig=as.numeric(apply(x$`Exposure:Station`[,2:3],1,prod)>=0)+1
  op=par(mar=c(4.2,9,3.8,2))
  plot(x,col=psig,yaxt="n")
  for (j in 1:length(psig)){
  axis(2,at=j,labels=rownames(x$`Exposure:Station`)[length(psig)-j+1],
       las=1,cex.axis=.8,col.axis=psig[length(psig)-j+1])
  }
  par(op)
}

sizmod$plors[[1]]

sizmod$ltmod[[1]]
## Post:Edge/Out       Pre:Mid      Post:Mid  Pre:Edge/Out 
##           "a"           "b"           "a"           "b"
plofun(sizmod$phmod[[1]])

sizmod$plors[[2]]

sizmod$ltmod[[2]]
## Post:Edge/Out       Pre:Mid      Post:Mid  Pre:Edge/Out 
##           "a"           "b"           "c"           "b"
plofun(sizmod$phmod[[2]])

sizmod$plors[[3]]

sizmod$ltmod[[3]]
## Post:Edge/Out       Pre:Mid      Post:Mid  Pre:Edge/Out 
##           "a"           "b"           "a"           "b"
plofun(sizmod$phmod[[3]])

sizmod$plors[[4]]

sizmod$ltmod[[4]]
## Post:Edge/Out       Pre:Mid      Post:Mid  Pre:Edge/Out 
##           "a"          "ab"           "c"           "b"
plofun(sizmod$phmod[[4]])

Alternative plots, these show the post-hoc comparisons of the locations from an ANOVA of treatment x location, see the above plots for all comparisons:

# Olympia
toplo <- sizmod %>% 
  mutate(
    phmod = purrr::map(phmod, function(x){
      x$`Exposure:Station` %>% 
        data.frame %>% 
        mutate(
          comp = row.names(.)
        )
    })
  ) %>% 
  select(Organism, var, phmod) %>% 
  unnest('phmod') %>% 
  filter(comp %in% c('Post:Edge/Out-Pre:Edge/Out', 'Post:Mid-Pre:Mid')) %>% 
  mutate(
    comp = factor(comp, 
      levels = c('Post:Edge/Out-Pre:Edge/Out', 'Post:Mid-Pre:Mid'), 
      labels = c('Out/Edge', 'Mid')
    ), 
    pcol = ifelse(p.adj <= 0.05, 'p < 0.05', 'ns'), 
    pcol = factor(pcol, levels = c('p < 0.05', 'ns')),
    var = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Delta~Length (mm)', 'Delta~Weight (g)')), 
    Organism = factor(Organism, levels = c('Olympia', 'Pacific'))
  ) 

p1 <- ggplot(toplo, aes(y = comp, x = diff, linetype = pcol, col = Organism)) + 
  geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
  geom_point(position = position_dodge(width = 0.3), size = 3) + 
  geom_errorbar(aes(xmin = lwr, xmax = upr), position = position_dodge(width = 0.3), width = 0.2) + 
  facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) + 
  scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
  theme_ipsum(plot_margin = margin(5, 5, 5, 5)) + 
  theme(
    strip.placement = 'outside',
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    strip.background = element_blank(), 
    axis.title.y = element_blank(), 
    legend.title = element_blank(), 
    legend.position = 'top',
    axis.title.x = element_text(hjust = 0.5),
    strip.text.x = element_text(hjust = 0.5),
    strip.text = element_text(size = 9),
    panel.background = element_rect()
    ) +
  labs(
    y = 'Location', 
    x = 'Effect size', 
    subtitle = '(a) Estimated means'
  )

# Olympia
toplo <- sizmod %>% 
  mutate(
    phmod = purrr::map(phmod, function(x){
      x$`Exposure:Station` %>% 
        data.frame %>% 
        mutate(
          comp = row.names(.)
        )
    })
  ) %>% 
  select(Organism, var, phmod) %>% 
  unnest('phmod') %>% 
  mutate(
    pcol = ifelse(p.adj <= 0.05, 'p < 0.05', 'ns'), 
    pcol = factor(pcol, levels = c('p < 0.05', 'ns')),
    var = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Delta~Length (mm)', 'Delta~Weight (g)')), 
    Organism = factor(Organism, levels = c('Olympia', 'Pacific'))
  ) 

p2 <- ggplot(toplo, aes(y = comp, x = diff, linetype = pcol, col = Organism)) + 
  geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
  geom_point(position = position_dodge(width = 0.3), size = 3) + 
  geom_errorbar(aes(xmin = lwr, xmax = upr), position = position_dodge(width = 0.3), width = 0.2) + 
  facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) + 
  scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
  theme_ipsum(plot_margin = margin(5, 5, 5, 5)) + 
  theme(
    strip.placement = 'outside',
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    strip.background = element_blank(), 
    axis.title.y = element_blank(), 
    legend.title = element_blank(), 
    legend.position = 'top',
    axis.title.x = element_text(hjust = 0.5),
    strip.text.x = element_text(hjust = 0.5),
    strip.text = element_text(size = 9), 
    panel.background = element_rect()
    ) +
  labs(
    y = 'Location', 
    x = 'Effect size', 
    subtitle = '(b) Multiple comparisons'
  )

p <- p1 + p2 + plot_layout(ncol = 1, heights = c(0.8, 1))

# pout
jpeg(here('figs/fieldcomp.jpeg'), height = 8, width = 7, family = 'serif', units = 'in', res = 500)
print(p)
dev.off()
knitr::include_graphics(here('figs/fieldcomp.jpeg'))

# size, species models
sizsppmod <- sizdat %>%
  mutate(
    Station = case_when(
      Station %in% c('EDG', 'OUT') ~ 'Edge/Out', 
      T ~ 'Mid'
    )
  ) %>% 
  group_by(Organism, var) %>% 
  mutate(val2 = as.numeric(scale(val))) %>% 
  pivot_longer(cols = val:val2, names_to = 'scl', values_to = 'val') %>% 
  group_by(var, scl) %>% 
  nest() %>% 
  mutate(
    lmmod = map(data, ~aov(val ~ Exposure * Station * Organism, data = .x)),
    phmod = map(lmmod, TukeyHSD),
    ltmod = map(phmod, function(x){
      modhsd <- x %>%
        .$`Exposure:Station:Organism` %>%
        data.frame
      pvals <- modhsd$p.adj
      names(pvals) <- rownames(modhsd)
      lets <- multcompLetters(pvals)

      return(lets)

    })
  )

toplo <- sizsppmod %>% 
  mutate(
    phmod = purrr::map(phmod, function(x){
      x$`Exposure:Station:Organism` %>% 
        data.frame %>% 
        mutate(
          comp = row.names(.)
        )
    })
  ) %>% 
  select(var, scl, phmod) %>% 
  unnest('phmod') %>% 
  mutate(
    pcol = ifelse(p.adj <= 0.05, 'p < 0.05', 'ns'), 
    pcol = factor(pcol, levels = c('p < 0.05', 'ns'))
  )

thm <- theme_ipsum() + 
  theme(
    strip.placement = 'outside',
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    strip.background = element_blank(), 
    axis.title.y = element_blank(), 
    legend.title = element_blank(), 
    legend.position = 'top',
    axis.title.x = element_text(hjust = 0.5),
    strip.text.x = element_text(hjust = 0.5),
    strip.text = element_text(size = 9), 
    panel.background = element_rect()
  )

toplo1 <- toplo %>% 
  filter(scl == 'val') %>% 
  mutate(
    var = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Delta~Length (mm)', 'Delta~Weight (g)'))
  )

p1 <- ggplot(toplo1, aes(y = comp, x = diff, linetype = pcol)) + 
  geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
  geom_point(size = 3) + 
  geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.2) + 
  facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) + 
  scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
  thm +
  labs(
    y = 'Location', 
    x = 'Effect size'
  )

toplo2 <- toplo %>% 
  filter(scl == 'val2') %>% 
  mutate(
    var = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Delta~Length (scaled)', 'Delta~Weight (scaled)')),
  )

p2 <- ggplot(toplo2, aes(y = comp, x = diff, linetype = pcol)) + 
  geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
  geom_point(size = 3) + 
  geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.2) + 
  facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) + 
  scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
  thm +
  labs(
    y = 'Location', 
    x = 'Effect size'
  )

# filter scaled post-hoc comparisons to those that are relevant
toplo3 <- toplo2 %>% 
  filter(comp %in% c("Pre:Mid:Pacific-Pre:Mid:Olympia", "Pre:Mid:Pacific-Post:Edge/Out:Pacific", 
                     "Post:Mid:Pacific-Pre:Mid:Pacific", "Post:Mid:Pacific-Pre:Edge/Out:Pacific", 
                     "Post:Mid:Pacific-Post:Mid:Olympia", "Post:Mid:Pacific-Post:Edge/Out:Olympia", 
                     "Post:Edge/Out:Pacific-Post:Mid:Olympia", "Post:Edge/Out:Pacific-Post:Edge/Out:Olympia")
           )

p3 <- ggplot(toplo3, aes(y = comp, x = diff, linetype = pcol)) + 
  geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
  geom_point(size = 3) + 
  geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.2) + 
  facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) + 
  scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
  thm +
  labs(
    y = 'Location', 
    x = 'Effect size'
  )

# pout
jpeg(here('figs/fieldcomp2.jpeg'), height = 7, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p1)
dev.off()

# pout
jpeg(here('figs/fieldcomp3.jpeg'), height = 7, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p2)
dev.off()

# pout
jpeg(here('figs/fieldcomp4.jpeg'), height = 5, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p3)
dev.off()
knitr::include_graphics(here('figs/fieldcomp2.jpeg'))

knitr::include_graphics(here('figs/fieldcomp3.jpeg'))

knitr::include_graphics(here('figs/fieldcomp4.jpeg'))

sclex <- sizsppmod %>% 
  select(-lmmod, -phmod, -ltmod) %>% 
  unnest('data') %>% 
  mutate(
    var1 = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Length (mm)', 'Weight (g)')),
    var2 = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Length (scaled)', 'Weight (scaled)'))
  ) 

thm <- theme_ipsum() + 
  theme(
    strip.placement = 'outside',
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    strip.background = element_blank(), 
    legend.title = element_blank(), 
    legend.position = 'top',
    axis.title.x = element_text(hjust = 0.5),
    axis.title.y = element_text(hjust = 0.5, size = 12),
    strip.text.x = element_text(hjust = 0.5),
    strip.text.y = element_text(hjust = 0.5),
    panel.spacing = unit(0.5, 'lines')
  )

toplo1 <- sclex %>% 
  filter(scl == 'val')

p1 <- ggplot(toplo1, aes(x = val, fill = Organism)) +
  geom_density(alpha = 0.7) + 
  facet_grid(Exposure ~ var1, switch = 'x', scales = 'free_x') + 
  thm + 
  labs(
    x = NULL, 
    y = 'Density', 
    subtitle = 'Measured values'
  )

toplo2 <- sclex %>% 
  filter(scl == 'val2')

p2 <- ggplot(toplo2, aes(x = val, fill = Organism)) +
  geom_density(alpha = 0.7) + 
  facet_grid(Exposure ~ var2, switch = 'x', scales = 'free_x') + 
  thm + 
  labs(
    x = NULL, 
    y = 'Density', 
    subtitle = 'Scaled values'
  )

# pout
jpeg(here('figs/fieldcomp5.jpeg'), height = 7, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p1)
dev.off()

jpeg(here('figs/fieldcomp6.jpeg'), height = 7, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p2)
dev.off()
knitr::include_graphics(here('figs/fieldcomp5.jpeg'))

knitr::include_graphics(here('figs/fieldcomp6.jpeg'))

Dissolution

p1 <- dissmod$plomod[[1]]
p2 <- dissmod$plomod[[2]]
p <- p2[[1]] + p2[[2]] + p1[[1]] + p1[[2]] + plot_layout(ncol = 2, guides = 'collect') & 
  theme(
    legend.position = 'bottom',
    plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), 'in')
    )

jpeg(here('figs/fielddiss.jpeg'), height = 7, width = 7, family = 'serif', units = 'in', res = 500)
print(p)
dev.off()
knitr::include_graphics(here('figs/fielddiss.jpeg'))

Pacific ANOVA:

dissmod %>% filter(Organism == 'Pacific') %>% pull(anomod) %>% .[[1]]
## Analysis of Variance Table
## 
## Response: Average.Score
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Exposure   2 13.481  6.7406  13.604 3.295e-05 ***
## Residuals 39 19.323  0.4955                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pacific contrasts:

dissmod %>% filter(Organism == 'Pacific') %>% pull(lmmod) %>% .[[1]] %>% estimate_contrasts('Exposure')
## Marginal Contrasts Analysis
## 
## Level1 | Level2 | Difference |        95% CI |   SE | t(39) |      p
## --------------------------------------------------------------------
## Edge   |    Mid |       0.99 | [ 0.29, 1.69] | 0.28 |  3.54 | 0.002 
## Out    |   Edge |       0.24 | [-0.48, 0.96] | 0.29 |  0.85 | 0.402 
## Out    |    Mid |       1.24 | [ 0.61, 1.86] | 0.25 |  4.96 | < .001
## 
## Marginal contrasts estimated at Exposure
## p-value adjustment method: Holm (1979)

Olympia ANOVA:

dissmod %>% filter(Organism == 'Olympia') %>% pull(anomod) %>% .[[1]]
## Analysis of Variance Table
## 
## Response: Average.Score
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Exposure   2 10.439  5.2194  11.362 0.0001423 ***
## Residuals 37 16.997  0.4594                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Olympia contrasts:

dissmod %>% filter(Organism == 'Olympia') %>% pull(lmmod) %>% .[[1]] %>% estimate_contrasts('Exposure')
## Marginal Contrasts Analysis
## 
## Level1 | Level2 | Difference |        95% CI |   SE | t(37) |      p
## --------------------------------------------------------------------
## Edge   |    Mid |       1.25 | [ 0.56, 1.95] | 0.28 |  4.52 | < .001
## Out    |   Edge |      -0.42 | [-1.11, 0.28] | 0.28 | -1.51 | 0.141 
## Out    |    Mid |       0.83 | [ 0.21, 1.46] | 0.25 |  3.37 | 0.004 
## 
## Marginal contrasts estimated at Exposure
## p-value adjustment method: Holm (1979)